#include "RBF.h"

using namespace std;

RBF::RBF(){loadConfig("default.txt");}

RBF::RBF(RBF* RB){
	this->hiddens = RB->hiddens;
	this->outputs = RB->outputs;
}

RBF::RBF(string file){
	loadConfig(file);
}

/**
	caso 0, se considera la matriz de covarianzas como la identidad
*/
void RBF::desvCase0( vector<vector<double> >& data){
	for(int i=0;i<this->nHiddens;i++)
		this->hiddens[i].desv_0 = 1;
}

/**
	caso1, matriz de covarianzas diagonal, multiplo de la identidad
	la varianza es el promedio de las varianzas en todas las direcciones
*/
void RBF::desvCase1( vector<vector<double> >& data){
	for(int i=0;i<this->nHiddens;i++){
		vector<double> desv;	
		int N = this->hide2data[i].size();	
		desv.resize(this->nInputs);
		for(int j=0;j<this->nInputs;j++){
			desv[j] = 0;
			for(int k=0; k<N;k++){
				desv[j] += pow( (data[hide2data[i][k]][j]-this->hiddens[i].mean[j]) , 2);
			}
			desv[j] /= (N-1);
			desv[j] = sqrt(desv[j]);
		}
		//view_vd(desv);
		this->hiddens[i].desv_1 = media(desv);
		cout<<"desv: "<<this->hiddens[i].desv_1<<endl;	
		cin.get();
	}
}

/**
	caso 2 (matriz de covarianza diagonal, con elementos distintos)
	cada matriz de covarianza tiene en la diagonal la respectiva varianza segun la direccion
*/
void RBF::desvCase2( vector<vector<double> >& data){
	for(int i=0;i<this->nHiddens;i++){
		int N = this->hide2data[i].size();	
		this->hiddens[i].desv_2.resize(this->nInputs);
		for(int j=0;j<this->nInputs;j++){
			this->hiddens[i].desv_2[j] = 0;
			for(int k=0; k<N;k++){
				this->hiddens[i].desv_2[j] += pow( (data[hide2data[i][k]][j]-this->hiddens[i].mean[j]) , 2);			
			}
			this->hiddens[i].desv_2[j] /= (N-1);	
			this->hiddens[i].desv_2[j] = sqrt(this->hiddens[i].desv_2[j]);
		}
		cout<<"desvio 2: "<<endl;
		view_vd(this->hiddens[i].desv_2);
		cin.get();
	}

}

/**
	caso 3, matriz de covarianza completa.
*/
void RBF::desvCase3( vector<vector<double> >& data){
	cout<<"en construccion desvCase3"<<endl;
	this->desvCase0(data);
}


/**
	caso simplificado 3 de las ppt, matriz identidad
*/
double RBF::calcOut0(int nGaussean, vector<double>& data){
	/*
	cout<<"dato: ";
	view_vd(data);
	cout<<"media: ";
	view_vd(this->hiddens[nGaussean].mean);
	*/
	return exp(-0.5*difference(data,this->hiddens[nGaussean].mean,2));
}

/**
	caso simplificado 2 de las ppt, matriz diagonal multiplo de la identidad
*/
double RBF::calcOut1(int nGaussean,vector<double>& data){
	double pi = 3.14159;
	double sigma = this->hiddens[nGaussean].desv_1;
	int N = this->nInputs;
	double exponent = -(1/(2*sigma*sigma))*difference(data,this->hiddens[nGaussean].mean,2);
	//double denominator = pow(2*pi,N) * sqrt(N) * sigma;
	//return (1/denominator)*exp(exponent);
	return exp(exponent);
}

/**
	caso simplificado 1 de las ppt, matriz diagonal general
*/
double RBF::calcOut2(int nGaussean,vector<double>& data){
	double pi = 3.14159;
	vector<double> sigma = this->hiddens[nGaussean].desv_2;
	vector<double> mu = this->hiddens[nGaussean].mean;
	int N = this->nInputs;
	double e = 0;
	for(int i=0; i<N; i++){
		e += (data[i]-mu[i])*(data[i]-mu[i]) / (sigma[i]*sigma[i]);
	}
	double exponent = -0.5*e;
	//double denominator = pow(2*pi,N) * norma_vec(sigma);
	//return (1/denominator)*exp(exponent);
	return exp(exponent);
}

/**
	caso matriz de covarianza completa.
*/
double RBF::calcOut3(int nGaussean,vector<double>& data){
	cout<<"en construccion calcOut3"<<endl;
	return this->calcOut0(nGaussean, data);
}


float RBF::uniform(){
	return ((double)rand() / ((double)(RAND_MAX)+(double)(1)));
}

void RBF::randomWeights(){
	for(int i=0;i<nOutputs;i++){
		for(int j=0;j<(nHiddens+1);j++){
			this->outputs[i].w[j] = uniform()-0.5;
		}	
	}
}

/**
	funcion wrapper para calcular la matriz de covarianza segun el tipo
*/
void RBF::calculateDesv( vector<vector<double> >& data){
	switch(this->typeDesv){
		case 3: this->desvCase3(data); break;
		case 2: this->desvCase2(data); break;
		case 1: this->desvCase1(data); break;
		default: this->desvCase0(data); break;
	}		
}

/**
	calcula la salida de una gausseana
	dependiendo el tipo de desvio utilizado se llama al correspondiente metodo
	esto es para no malgastar recursos generalizando
*/
double RBF::gausseanOut(int nGaussean, vector<double>& data){
	double ret;
	switch(this->typeDesv){
		case 3: ret = this->calcOut3(nGaussean, data); break;
		case 2: ret = this->calcOut2(nGaussean, data); break;
		case 1: ret = this->calcOut1(nGaussean, data); break;
		default: ret = this->calcOut0(nGaussean, data); break;
	}
	return ret;
}

//------------------------------------------------------------
void RBF::batchKmeans(int nData, vector<vector<double> >& data){
	cout<<"En batch..."<<endl;
	//Asigna cada patrón de entrada a un conjunto aleatoriamente
	for(int i=0; i<nData; i++){
		int j = int(rand() % this->nHiddens);	//Genera el índice del conjunto aleatoriamente
		this->hide2data[j].push_back(i);	//Ingresa el índice dato al conjunto
		this->data2hide[i] = j;
	}
	//view_vi(this->data2hide);
	if(this->plot)
		this->plotSets(data);
	//cout<<"asigna centroides random."<<endl;
	bool next = true;
	while(next){
		//Calculo y traslado de centroides
		for(int i=0; i<nHiddens; i++){	//Recorre conjuntos
			vector<double> means;
			means.resize(nInputs);
			for(int j=0; j<nInputs; j++){	//Recorre dimensiones de entrada
				int h2dLength = this->hide2data[i].size();	//Número de datos en el conjunto
				means[j]=0;
				for(int k=0; k<h2dLength; k++){	//Recorre patrones por conjunto en cada dimensión
					means[j]+=data[this->hide2data[i][k]][j];//Suma los patrones
				}
				if(h2dLength!=0)
					means[j]/=h2dLength;	//Promedio
			}
			this->hiddens[i].mean = means;	//Reasigna los centroides
			//cout<<"media "<<i<<endl;
			view_vd(this->hiddens[i].mean);
		}
		next = false;
		//Reasignación de patrones a conjuntos
		for(int i=0; i<nData; i++){	//Recorre datos de entrada
			int indexMin = 0;	//Indice conjunto más cercano
			double min = INT_MAX;	//Distancia mínima
			for(int j=0; j<this->nHiddens; j++){	//Recorre filas de conjuntos
				double d = euclideanDist(nInputs,data[i],this->hiddens[j].mean);//Distancia entre patron y centroides
				if(d < min){
					indexMin = j;
					min = d;
				}
			}
			if(indexMin != data2hide[i]){	//Si el patron cambia de centroide
				next = true;
				vector<int>::iterator it;
				it = find(this->hide2data[this->data2hide[i]].begin(), this->hide2data[this->data2hide[i]].end(),i);	//Encuentra el indice del patron en el conjunto viejo
				this->hide2data[this->data2hide[i]].erase(it);	//Lo borra
				this->hide2data[indexMin].push_back(i);	//Lo agrega al conjunto al que pertenece ahora
				this->data2hide[i] = indexMin;	//Actualiza el conjunto para su respectivo patrón
			}
		}
		if(this->plot)
			this->plotSets(data);
	}
	cout<<"termina batch ..."<<endl;
}
//------------------------------------------------------------
void RBF::onlineKmeans(int nData, vector<vector<double> >& data){
	cout<<"En online k-means"<<endl;
	for(int i=0; i<this->nHiddens; i++){
		int index;
		index = int(rand() % nData);
		this->hiddens[i].mean = data[index];//Corroborar que el índice no se repita
	}
	int epoca = 1;
	bool next = true;
	double JBefore = 0;
	while(next){
		if(this->plot){ //solo si voy a plotear
			for(int i=0; i<this->nHiddens;i++){
				this->hide2data[i].clear();
			}
		}
		double J = 0; //contendra la suma de las distancias de todos los patrones a su media
		//para todos los patrones encuentro la media mas cercana a c/u y la actualizo 
		for(int i=0; i<nData; i++){
			int indexMin = 0;	//Indice conjunto más cercano
			double min = INT_MAX;	//Distancia mínima
			for(int j=0; j<this->nHiddens; j++){	//Recorre filas de conjuntos
				double d = euclideanDist(nInputs,data[i],this->hiddens[j].mean);	//Distancia entre patron y centroides
				if(d < min){
					indexMin = j;
					min = d;
				}
			}
			for(int j=0; j<nInputs; j++){
				float delta = etaOnline*(data[i][j]-this->hiddens[indexMin].mean[j]);
				this->hiddens[indexMin].mean[j] += delta;
			}
		}

		//calculo el error luego de todas las actualizaciones
		double distance;
		//double error = 0;
		for(int i=0; i<nData; i++){
			int indexMin = 0;	//Indice conjunto más cercano
			double min = INT_MAX;	//Distancia mínima
			for(int j=0; j<this->nHiddens; j++){	//Recorre filas de conjuntos
				distance = euclideanDist(nInputs,data[i],this->hiddens[j].mean);	//Distancia entre patron y centroides
				if(distance < min){
					indexMin = j;
					min = distance;
				}
			}
			if(this->plot){ //aca me sirve solo si ploteo
				this->hide2data[indexMin].push_back(i);	//Ingresa el índice dato al conjunto
		  	this->data2hide[i] = indexMin;
			}
			J += distance;
		}
		if(fabs(J-JBefore)<deltaError){
			next = false;
		}
		if(epoca > 1000){ //parametrizar max epocas
			next = false;
		}
		JBefore = J;
		
		if(this->plot)
			this->plotSets(data);

		epoca++;
	}

	//luego de entrenar tengo que asignar patrones a gausseanas para calcular desvios
	// en batch esto se va asignando durante el propio entrenamiento
	for(int i=0; i<this->nHiddens;i++){
		this->hide2data[i].clear();
	}
	for(int i=0; i<nData; i++){
		int indexMin = 0;	//Indice conjunto más cercano
		double min = INT_MAX;	//Distancia mínima
		for(int j=0; j<this->nHiddens; j++){	//Recorre filas de conjuntos
			double d = euclideanDist(nInputs,data[i],this->hiddens[j].mean);	//Distancia entre patron y centroides
			if(d < min){
				indexMin = j;
				min = d;
			}
		}
		this->hide2data[indexMin].push_back(i);	//Ingresa el índice dato al conjunto
		this->data2hide[i] = indexMin;
	}

	cout<<"Salida de online Kmeans en la iteracion numero: "<<epoca<<endl<<endl;

}
//------------------------------------------------------------
void RBF::lms(int nData, vector<vector<double> >& data, vector<vector<int> >& targets){
	cout<<"comienza lms ..."<<endl;
	double delta, error, sum, E, Eav, EavSigno;
	int ESigno;
	bool errorSigno;
	this->randomWeights();
	int epoca = 1;
	bool next = true;

	//solo requiero calcular una vez las salidas de las capas ocultas (estan estaticas)
	for(int j=0;j<this->nHiddens;j++){
		this->hiddens[j].out.resize(nData);
	}
	for(int i=0;i<nData;i++){
 		//calculo la salida en la capa oculta
		for(int j=0;j<this->nHiddens;j++){
			this->hiddens[j].out[i] = this->gausseanOut(j,data[i]);	
		}
	}
	while(next){
		Eav = 0; //error promedio para la epoca
		EavSigno = 0;
		for(int i=0;i<nData;i++){			//cantidad de patrones
			//calculo la salida en la capa de salida y corrijo los pesos
			for(int j=0;j<this->nOutputs;j++){	//cantidad de neuronas de salida
				sum = (-1)*this->outputs[j].w[0]; //bias
				for(int k=0;k<this->nHiddens;k++){
					sum += 	this->outputs[j].w[k+1] * this->hiddens[k].out[i];//prod punto entre la salida de la oculta y los pesos
				}
				
				error = sum - targets[i][j]; //error para esta neurona
				for(int k=0;k<(this->nHiddens+1);k++){
					if(k>0)
						delta = this->etaLMS*error*this->hiddens[k-1].out[i];
					else
						delta = this->etaLMS*error*-1;	
					this->outputs[j].w[k] -= delta;
				}
			}
		}//end for i
		//luego de recorrer todos los patrones y actualizar, calculo el error para esta epoca
		ESigno = 0;
		for(int i=0;i<nData;i++){
			E = 0;
			errorSigno = false; //bandera
			//Calculo las salidas de la capa de salida
			for(int j=0;j<this->nOutputs;j++){
				sum = (-1)*this->outputs[j].w[0];
				for(int k=0;k<this->nHiddens;k++){
					sum += 	this->outputs[j].w[k+1] * this->hiddens[k].out[i];//prod punto entre la salida de la oculta y los pesos
				}
				if(sign(sum)!=targets[i][j]) //si hay error de clasificacion
					errorSigno = true;
				error = sum - targets[i][j];
				E += 0.5*error*error;
			}
			Eav += E;
			if(errorSigno){
				ESigno++;
			}
		}//end for i
		
		Eav /= nData;
	  cout<<"Num Errores de signo "<<ESigno<<" de "<<nData<<endl;
		if(ESigno>0)
			EavSigno = (double)ESigno/(double)nData;	
		else
			EavSigno = 0.0;
		cout<<"Numero de epoca: "<<epoca;
		cout<<"  - Error: "<<Eav<<endl;
		//criterios de parada
		if(EavSigno == 0){
			next = false;
			cout<<endl<<"-------------------------------------------------------"<<endl;
			cout<<"Salida por Clasificacion correcta"<<endl<<"Numero de epocas: "<<epoca<<endl;
			cout<<"Error de clasificacion pasando por función signo: "<<EavSigno*100<<" %"<<endl;
			cout<<"Error de clasificacion: "<<Eav<<endl;
		}
		if(Eav < this->deltaError){
			next = false;
			cout<<endl<<"-------------------------------------------------------"<<endl;
			cout<<"Salida por delta de error"<<endl<<"|Eav|: "<<Eav<<endl<<"Numero de epocas: "<<epoca<<endl;
			cout<<"Error de clasificacion pasando por función signo: "<<EavSigno*100<<" %"<<endl;
			cout<<"Error de clasificacion: "<<Eav<<endl;
		}
		if(epoca >= maxEpocas){
			next = false;
			cout<<endl<<"-------------------------------------------------------"<<endl;
			cout<<"Salida por maximo de epocas: "<<epoca<<endl;
			cout<<"Error de clasificacion pasando por función signo "<<EavSigno*100<<" %"<<endl;
			cout<<"Error de clasificacion: "<<Eav<<endl;
			cout<<endl<<"-------------------------------------------------------"<<endl;
		}
		epoca++;
	}//end while
}

//------------------------------------------------------------
void RBF::training(vector<vector<double> >& data, vector<vector<int> >& targets){
	srand(time(NULL));
	//cout<<"nHiddens: "<<nHiddens<<endl<<"nOutputs: "<<nOutputs<<endl<<"nInputs: "<<nInputs<<endl<<"typeDesv: "<<typeDesv<<endl<<"maxEpocas: "<<maxEpocas<<endl;
	int nData = data.size();
	this->data2hide.resize(nData);
	if(this->typeKmeans)
		this->batchKmeans((int)nData,data);
	else
		this->onlineKmeans((int)nData,data);
	this->calculateDesv(data);
	this->lms(nData,data,targets);		
}

//------------------------------------------------------------
double RBF::test(vector<vector<double> >& data, vector<vector<int> >& targets){
	int nData = data.size();
	double Eav, sum;
	int E = 0;
	bool error;
	for(int j=0;j<this->nHiddens;j++){
		this->hiddens[j].out.resize(nData);
	}
	//calculo la salida en capas intermedias
	for(int i=0;i<nData;i++){
		for(int j=0;j<this->nHiddens;j++){
			this->hiddens[j].out[i] = this->gausseanOut(j,data[i]);	
		}
	}

	E = 0; //acumulo los errores de clasificacion
	for(int i=0;i<nData;i++){
		error = false;
		for(int j=0;j<this->nOutputs;j++){
			sum = (-1)*this->outputs[j].w[0];
			for(int k=0;k<this->nHiddens;k++){
				sum += 	this->outputs[j].w[k+1] * this->hiddens[k].out[i];
			}
			if(sign(sum)!=targets[i][j])
				error = true;
		}
		//si hubo error de clasificacion, sumo
		if(error)
			E++;
	}//end for i
	cout<<"Cantidad de errores en test "<<E<<" de "<<nData<<endl;
	if(E>0)
		Eav = (double)E/(double)nData; //porcentaje de error	
	else
		Eav = 0.0;
	cout<<"Porcentaje de error "<<Eav*100<<" %"<<endl;
	return Eav;
}

//------------------------------------------------------------
void RBF::loadConfig(string file){
	ifstream archi;
	archi.open(file.c_str());

	if(archi.is_open()){
		string str_line;
		char* line;
		//Titulo
		getline(archi,str_line);
		//Número de entradas
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->nInputs = atoi(line);
		//Numero de neuronas ocultas
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->nHiddens = atoi(line);
		//Numero de neuronas de salida
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->nOutputs = atoi(line);
		//Tipo de K means
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->typeKmeans = atoi(line);
		//Tipo de desvío
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->typeDesv = atoi(line);
		//Tasa de aprendizaje online K means
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->etaOnline = atof(line);
		//Tasa de aprendizaje LMS
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->etaLMS = atof(line);
		//Delta error de entrenamiento
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->deltaError = atof(line);
		//maximo de epocas de entrenamiento en supervisado
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->maxEpocas = atoi(line);
		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->plot = atoi(line);
		}		
	archi.close();
	
	this->hiddens.resize(this->nHiddens);
	this->outputs.resize(this->nOutputs);
	this->hide2data.resize(nHiddens);	//Redimensiona la cantidad de filas (cantidad de conjuntos)

	for(int i=0; i<this->nHiddens; i++){
		this->hiddens[i].mean.resize(this->nInputs);
		switch(this->typeDesv){
			case 3: this->hiddens[i].desv_3.resize(this->nInputs); break; //falta resize para la 2da dimension, an no se hace xq no se implementa matriz llena
			case 2: this->hiddens[i].desv_2.resize(this->nInputs); break;
		}
	}

	for(int i=0; i<this->nOutputs; i++){
		this->outputs[i].w.resize(this->nHiddens+1);
	}
}


//ploteo los conjuntos tal como van quedando
		//solo si estamos en 2 dimensiones en el espacio de entradas		
void RBF::plotSets(vector<vector<double> >& data){
	if(nInputs == 2){		
		vector<string> paths;
		for(unsigned int i=0;i<this->hide2data.size();i++){
			vector<vector<double> > aux;
			for(unsigned int j=0; j<this->hide2data[i].size();j++){
				aux.push_back(data[this->hide2data[i][j]]);			
			}	
			stringstream cvz;
			cvz<<i;
			string name = "Plots/rbf_gaussean_" + cvz.str() + ".txt";
			if(aux.size()>0){
				saveFileX(aux,name);
				paths.push_back(name);
			}
		}
		this->printData(paths);
	}
}

//-------------------------
//Global para poner replotear en la misma instancia!
GNUplot gp; 
//------------------------

/**
Printeamos las clases con GNUPLOT, pasando a traves de un arreglo
*/
void RBF::printData(vector<string> paths) { 
	stringstream ss;
	ss.flags(std::ios::fixed);
	ss.precision(4);
	ss.width(0);
	int lt[]={-1,1,2,3,4};
	ss<<"plot [-13:13][-13:13]";
	for (unsigned int i=0;i<paths.size();i++){
		ss<<" '"<<paths[i].c_str()<<"' lt "<<lt[i%5];
		if(i+1 < paths.size())
			ss<<",";
		else
			ss<<" ;";
	}
	//cout<<ss.str()<<endl;
	gp(ss.str());
	for (int i=0;i<100000000;i++){ /*algo*/}
}
